********************************************************************************
* Estimate HSV on total income
********************************************************************************

include choose_version.do

putexcel set $figure_path/cps_results_`version_robust'.xlsx, sheet(tax_functions) modify open
putexcel A1 = "Estimation"
putexcel B1 = "lda (HSV)"
putexcel C1 = "tau (HSV)"
putexcel D1 = "lda (FGNV)"
putexcel E1 = "theta (FGNV)"
putexcel F1 = "m (FGNV)"
putexcel G1 = "xi (FGNV)"
putexcel H1 = "T"
putexcel I1 = "Deviation level"
putexcel J1 = "Deviation log"

* Generate additional income variables
gen inc_tot_hh_at = inc_tot_hh - tax_transfer
gen log_inc_tot_hh = log(inc_tot_hh)
gen log_inc_tot_hh_at = log(inc_tot_hh_at)

* Mean pre-tax total income (to be used as normalization for all tax-transfer function estimates)
qui sum inc_tot_hh [fw = asecwth], d
scalar define ybar = r(mean)
	
* Linear estimate in logs
reg log_inc_tot_hh_at log_inc_tot_hh [fw = asecwth]						// linear regression
predict double log_inc_tot_hh_at_pred_hsv_log							// predicted log after tax income of HSV estimate in logs
gen inc_tot_hh_at_pred_hsv_log = exp(log_inc_tot_hh_at_pred_hsv_log)	// predicted after tax income of HSV estimate in logs
	
* Coeffients
display "tau = " 1-_b[log_inc_tot_hh]
display "lda_hsv = " exp(_b[_cons])*ybar^(-(1-_b[log_inc_tot_hh]))
putexcel A2 = "HSV OLS Logs"
putexcel B2 = (exp(_b[_cons])*ybar^(-(1-_b[log_inc_tot_hh])))
putexcel C2 = 1-_b[log_inc_tot_hh] 
putexcel D2 = "-"
putexcel E2 = "-"
putexcel F2 = "-"
putexcel G2 = "-"
putexcel H2 = "-"

* Deviations
gen hsv_log_ssd_dollar = (inc_tot_hh_at-inc_tot_hh_at_pred_hsv_log)^2
sum hsv_log_ssd_dollar [fw = asecwth]
putexcel I2 = (sqrt(`r(mean)'))
gen hsv_log_ssd_log = (log(inc_tot_hh_at)-log(inc_tot_hh_at_pred_hsv_log))^2
sum hsv_log_ssd_log [fw = asecwth]
putexcel J2 = (sqrt(`r(mean)'))

* Nonlinear estimate in logs
nl (log_inc_tot_hh_at = log({lda_hsv=0.8}) + (1-{tau=0.2})*log_inc_tot_hh + {tau=0.2}*log(ybar)) [fw = asecwth]
*nl (log_inc_tot_hh_at = log({lda_hsv=0.8}) + (1-{tau=0.2})*log_inc_tot_hh + {tau=0.2}*log(ybar)) [fw = asecwth], hasconstant(lda_hsv)	// same result as previous
*nl (log_inc_tot_hh_at = log({lda_hsv=0.8}) + (1-{tau=0.2})*log_inc_tot_hh + {tau=0.2}*log(ybar)) [fw = asecwth], nocons // same result as previous (though it should have a constant)
putexcel A3 = "HSV NL Logs"
putexcel B3 = (r(table)[1,1])
putexcel C3 = (r(table)[1,2])
putexcel D3 = "-"
putexcel E3 = "-"
putexcel F3 = "-"
putexcel G3 = "-"
putexcel H3 = "-"

* Nonlinear estimate in levels
nl (inc_tot_hh_at = {lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2})*inc_tot_hh) [fw = asecwth]
*nl (inc_tot_hh_at = {lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2})*inc_tot_hh) [fw = asecwth], nocons	// same result as previous
predict double inc_tot_hh_at_pred_hsv_level		// predicted after tax income of HSV estimate in levels
putexcel A4 = "HSV NL Levels"
putexcel B4 = (r(table)[1,1])
putexcel C4 = (r(table)[1,2])
putexcel D4 = "-"
putexcel E4 = "-"
putexcel F4 = "-"
putexcel G4 = "-"
putexcel H4 = "-"

* Deviations
gen hsv_lev_ssd_dollar = (inc_tot_hh_at-inc_tot_hh_at_pred_hsv_level)^2
sum hsv_lev_ssd_dollar [fw = asecwth]
putexcel I4 = (sqrt(`r(mean)'))
gen hsv_lev_ssd_log = (log(inc_tot_hh_at)-log(inc_tot_hh_at_pred_hsv_level))^2
sum hsv_lev_ssd_log [fw = asecwth]
putexcel J4 = (sqrt(`r(mean)'))

* Nonlinear estimate in levels: matching tax-transfer payment instead of 
nl (tax_transfer = inc_tot_hh*(1-{lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2}))) [fw = asecwth]				// same coefficients as previous
*nl (tax_transfer = inc_tot_hh*(1-{lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2}))) [fw = asecwth], nocons		// same result as previous
predict double tax_transfer_pred_hsv_level		// predicted tax-transfer payment (receipt) of HSV estimate in levels

********************************************************************************
* Estimate HSV plus intercept on total income
********************************************************************************

* Nonlinear estimate in levels
nl (inc_tot_hh_at = {lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2})*inc_tot_hh + {T=3000}) [fw = asecwth] 					// start from something other than 0
*nl (inc_tot_hh_at = {lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2})*inc_tot_hh + {T=3000}) [fw = asecwth], hasconstant(T) 	// same as previous
predict double inc_tot_hh_at_pred_hsvT_level	// predicted after tax income of HSV+T estimate in levels					
putexcel A5 = "HSV+T NL Levels"
putexcel B5 = (r(table)[1,1])
putexcel C5 = (r(table)[1,2])
putexcel D5 = "-"
putexcel E5 = "-"
putexcel F5 = "-"
putexcel G5 = "-"
putexcel H5 = (r(table)[1,3])

* Deviations
gen hsvT_lev_ssd_dollar = (inc_tot_hh_at-inc_tot_hh_at_pred_hsvT_level)^2
sum hsvT_lev_ssd_dollar [fw = asecwth]
putexcel I5 = (sqrt(`r(mean)'))
gen hsvT_lev_ssd_log = (log(inc_tot_hh_at)-log(inc_tot_hh_at_pred_hsvT_level))^2
sum hsvT_lev_ssd_log [fw = asecwth]
putexcel J5 = (sqrt(`r(mean)'))

* Nonlinear estimate in logs
nl (log_inc_tot_hh_at = log({lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2})*inc_tot_hh + {T=3000})) [fw = asecwth] 
*nl (log_inc_tot_hh_at = log({lda_hsv=0.8}*(inc_tot_hh/ybar)^(-{tau=0.2})*inc_tot_hh + {T=3000})) [fw = asecwth], nocons
putexcel A6 = "HSV+T NL Logs"
putexcel B6 = (r(table)[1,1])
putexcel C6 = (r(table)[1,2])
putexcel D6 = "-"
putexcel E6 = "-"
putexcel F6 = "-"
putexcel G6 = "-"
putexcel H6 = (r(table)[1,3])

********************************************************************************
* Estimate FGNV tax function on labor income
********************************************************************************

* Labor income minus taxes
gen inc_lab_hh_at = inc_lab_hh - tax

* Nonlinear estimate in levels on income
nl (inc_lab_hh_at = inc_lab_hh*(1-exp(log({lda_fgnv=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3})))) [fw = asecwth]
*nl (inc_lab_hh_at = inc_lab_hh*(1-exp(log({lda_fgnv=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3})))) [fw = asecwth], nocons	// same as previous
putexcel A7 = "FGNV Tax NL Levels"
putexcel B7 = "-"
putexcel C7 = "-"
putexcel D7 = (r(table)[1,1])
putexcel E7 = (r(table)[1,2])
putexcel F7 = "-"
putexcel G7 = "-"
putexcel H7 = "-"

* Nonlinear estimate in logs
gen log_inc_lab_hh_at = log(inc_lab_hh_at)
nl (log_inc_lab_hh_at = log(inc_lab_hh*(1-exp(log({lda_fgnv=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3}))))) [fw = asecwth]

* Nonlinear estimate in levels on tax payment
nl (tax = exp(log({lda_fgnv=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3}))*inc_lab_hh) [fw = asecwth]			// same coefficients as estimates on income in levels
*nl (tax = exp(log({lda_fgnv=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3}))*inc_lab_hh) [fw = asecwth], nocons	// same as previous
predict double tax_pred_fgnv_lev	// predicted tax of FGNV tax estimate in levels		

* Estimate on rates
nl (tax_rate = exp(log({lda_fgnv=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3}))) [fw = asecwth]
predict double tax_rate_fgnv		// predicted tax rate of FGNV tax rate estimate

********************************************************************************
* Estimate FGNV transfer function on total income
********************************************************************************

* Nonlinear estimate in levels on transfer payment
nl (transfer = {m=0.2}*ybar*(2*exp(-{xi=3}*(inc_tot_hh/ybar)))/(1+exp(-{xi=3}*(inc_tot_hh/ybar)))) [fw = asecwth]
*nl (transfer = {m=0.2}*ybar*(2*exp(-{xi=3}*(inc_tot_hh/ybar)))/(1+exp(-{xi=3}*(inc_tot_hh/ybar)))) [fw = asecwth], nocons	// same as previous
predict double transfer_pred_fgnv_lev	// predicted transfer of FGNV transfer estimate in levels		
putexcel A8 = "FGNV Transfer NL Levels"
putexcel B8 = "-"
putexcel C8 = "-"
putexcel D8 = "-"
putexcel E8 = "-"
putexcel F8 = (r(table)[1,1])
putexcel G8 = (r(table)[1,2])
putexcel H8 = "-"

* Nonlinear estimate in logs
gen log_transfer = log(transfer)
nl (log_transfer = log({m=0.2}*ybar*(2*exp(-{xi=3}*(inc_tot_hh/ybar)))/(1+exp(-{xi=3}*(inc_tot_hh/ybar))))) [fw = asecwth]

* Estimate on rates
nl (transfer_rate = ({m=0.2}*ybar*(2*exp(-{xi=3}*(inc_tot_hh/ybar)))/(1+exp(-{xi=3}*(inc_tot_hh/ybar)))/inc_tot_hh)) [fw = asecwth]

********************************************************************************
* Estimate FGNV tax-transfer function on total income
********************************************************************************

* Nonlinear estimate jointly in levels
nl (inc_tot_hh_at = inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3})) + {m=0.2}*ybar*(2*exp(-{xi=3}*inc_tot_hh/ybar))/(1+exp(-{xi=3}*inc_tot_hh/ybar))) [fw = asecwth]
*nl (inc_tot_hh_at = inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.3})) + {m=0.2}*ybar*(2*exp(-{xi=3}*inc_tot_hh/ybar))/(1+exp(-{xi=3}*inc_tot_hh/ybar))) [fw = asecwth], nocons	// same as previous
predict double inc_tot_hh_at_pred_fgnv_lev		// predicted after tax income of joint FGNV estimate in levels
putexcel A9 = "FGNV t&T NL Levels"
putexcel B9 = "-"
putexcel C9 = "-"
putexcel D9 = (r(table)[1,1])
putexcel E9 = (r(table)[1,2])
putexcel F9 = (r(table)[1,3])
putexcel G9 = (r(table)[1,4])
putexcel H9 = "-"

gen transfer_pred_fgnv_lev_joint = r(table)[1,3]*ybar*(2*exp(-r(table)[1,4]*inc_tot_hh/ybar))/(1+exp(-r(table)[1,4]*inc_tot_hh/ybar))	// predicted transfer of FGNV joint estimate in levels
gen tax_pred_fgnv_lev_joint = inc_lab_hh*exp(log(r(table)[1,1])*(inc_lab_hh/ybar)^(-2*r(table)[1,2]))									// predicted tax of FGNV joint estimate in levels

* Deviations
gen fgnv_joint_ssd_dollar = (inc_tot_hh_at-inc_tot_hh_at_pred_fgnv_lev)^2
sum fgnv_joint_ssd_dollar [fw = asecwth]
putexcel I9 = (sqrt(`r(mean)'))
gen fgnv_joint_ssd_log = (log(inc_tot_hh_at)-log(inc_tot_hh_at_pred_fgnv_lev))^2
sum fgnv_joint_ssd_log [fw = asecwth]
putexcel J9 = (sqrt(`r(mean)'))

* Nonlinear estimate jointly in logs
nl (log_inc_tot_hh_at = log(inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.2})) + {m=0.15}*ybar*(2*exp(-{xi=2}*inc_tot_hh/ybar))/(1+exp(-{xi=2}*inc_tot_hh/ybar)))) [fw = asecwth]
*nl (log_inc_tot_hh_at = log(inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.2})) + {m=0.15}*ybar*(2*exp(-{chi=2}*inc_tot_hh/ybar))/(1+exp(-{chi=2}*inc_tot_hh/ybar)))) [fw = asecwth], nocons // same as previous
predict double log_inc_tot_hh_at_pred_fgnv_log	// predicted log after tax income of joint FGNV estimate in logs 
putexcel A10 = "FGNV t&T NL Logs"
putexcel B10 = "-"
putexcel C10 = "-"
putexcel D10 = (r(table)[1,1])
putexcel E10 = (r(table)[1,2])
putexcel F10 = (r(table)[1,3])
putexcel G10 = (r(table)[1,4])
putexcel H10 = "-"

********************************************************************************
* Estimate FGNV function with transfer restricted to not phase out
********************************************************************************

* Restrict xi = 0
nl (inc_tot_hh_at = inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.2})) + {m=0.15}*ybar) [fw = asecwth]
*nl (inc_tot_hh_at = inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.2})) + {m=0.15}*ybar) [fw = asecwth], nocons	// same as previous
putexcel A11 = "FGNV+m NL Levels"
putexcel B11 = "-"
putexcel C11 = "-"
putexcel D11 = (r(table)[1,1])
putexcel E11 = (r(table)[1,2])
putexcel F11 = (r(table)[1,3])
putexcel G11 = "-"
putexcel H11 = "-"

* Estimate with T instead of m
nl (inc_tot_hh_at = inc_tot_hh - inc_lab_hh*exp(log({lda=0.2})*(inc_lab_hh/ybar)^(-2*{theta=0.2})) + {T=1000}) [fw = asecwth]
putexcel A12 = "FGNV+T NL Levels"
putexcel B12 = "-"
putexcel C12 = "-"
putexcel D12 = (r(table)[1,1])
putexcel E12 = (r(table)[1,2])
putexcel F12 = "-"
putexcel G12 = "-"
putexcel H12 = (r(table)[1,3])

putexcel save

********************************************************************************
* Save results/predictions of estimates
********************************************************************************

* Predicted after tax income of different variants
gen inc_tot_hh_at_pred_fgnv_lev_sep = inc_tot_hh + transfer_pred_fgnv_lev - tax_pred_fgnv_lev	// after tax income predicted using separate FGNV tax and transfer estimates
gen inc_tot_hh_at_pred_hsv_lev = inc_tot_hh - tax_transfer_pred_hsv_level						// after tax income predicted using HSV in levels estimate
gen inc_tot_hh_at_pred_fgnv_log = exp(log_inc_tot_hh_at_pred_fgnv_log)							// after tax income predicted using FGNV in logs estimate

* Deviations
gen fgnv_sep_ssd_dollar = (inc_tot_hh_at-inc_tot_hh_at_pred_fgnv_lev_sep)^2
sum fgnv_sep_ssd_dollar [fw = asecwth]
putexcel I7 = (sqrt(`r(mean)'))
putexcel I8 = (sqrt(`r(mean)'))
gen fgnv_sep_ssd_log = (log(inc_tot_hh_at)-log(inc_tot_hh_at_pred_fgnv_lev_sep))^2
sum fgnv_sep_ssd_log [fw = asecwth]
putexcel J7 = (sqrt(`r(mean)'))
putexcel J8 = (sqrt(`r(mean)'))

* In logs 
gen log_inc_tot_hh_at_pred_fgnv_lev = log(inc_tot_hh_at_pred_fgnv_lev) 							// log after tax income predicted using separate FGNV tax and transfer estimates
gen log_inc_tot_hh_at_pred_fgnv_levs = log(inc_tot_hh_at_pred_fgnv_lev_sep)						// log after tax income predicted using HSV in levels estimate
gen log_inc_tot_hh_at_pred_hsv_lev = log(inc_tot_hh_at_pred_hsv_lev)							// log after tax income predicted using FGNV in logs estimate

* Save results of estimates
xtile inc_tot_hh_groups = inc_tot_hh [fw = asecwth], nq(40)
xtile inc_lab_hh_groups = inc_lab_hh [fw = asecwth], nq(40)

putexcel set $figure_path/cps_results_`version_robust'.xlsx, sheet(tax_functions_figures) modify open
putexcel A1 = "Group"
putexcel B1 = "Pre-tax income"
putexcel C1 = "Post-tax income"
putexcel D1 = "Log pre-tax income"
putexcel E1 = "Log post-tax income"
putexcel F1 = "HSV logs: log post-tax income"
putexcel G1 = "HSV logs: post-tax income"
putexcel H1 = "HSV levels: post-tax income"
putexcel I1 = "HSV+T levels: post-tax income"
putexcel J1 = "FGNV levels (comb): post-tax income"
putexcel K1 = "FGNV levels (sep): post-tax income"
putexcel L1 = "FGNV tax"
putexcel M1 = "Tax"
putexcel N1 = "FGNV transfer"
putexcel O1 = "Transfer"
putexcel P1 = "FGNV transfer (joint est)"
putexcel Q1 = "FGNV tax (joint est)"
putexcel R1 = "SNAP"
putexcel S1 = "Housing Assistance"
putexcel T1 = "Welfare"
putexcel U1 = "Credits"
putexcel V1 = "State Credits"
putexcel W1 = "Pre-tax labor income"
putexcel X1 = "FGNV tax by lab. inc. group"
putexcel Y1 = "Tax by lab. inc. group"

local counter_aux = 1
forvalues i_group = 1(1)40 {
	
	local counter_aux = `counter_aux'+1
	putexcel A`counter_aux' = `i_group'
	
	sum inc_tot_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel B`counter_aux' = `r(mean)'
	
	sum inc_tot_hh_at if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel C`counter_aux' = `r(mean)'
	
	sum log_inc_tot_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel D`counter_aux' = `r(mean)'
	
	sum log_inc_tot_hh_at if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel E`counter_aux' = `r(mean)'
	
	sum log_inc_tot_hh_at_pred_hsv_log if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel F`counter_aux' = `r(mean)'
	
	sum inc_tot_hh_at_pred_hsv_log if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel G`counter_aux' = `r(mean)'
	
	sum inc_tot_hh_at_pred_hsv_level if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel H`counter_aux' = `r(mean)'
	
	sum inc_tot_hh_at_pred_hsvT_level if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel I`counter_aux' = `r(mean)'
	
	sum inc_tot_hh_at_pred_fgnv_lev if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel J`counter_aux' = `r(mean)'
	
	sum inc_tot_hh_at_pred_fgnv_lev_sep if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel K`counter_aux' = `r(mean)'
	
	sum tax_pred_fgnv_lev if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel L`counter_aux' = `r(mean)'
	
	sum tax if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel M`counter_aux' = `r(mean)'
	
	sum transfer_pred_fgnv_lev if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel N`counter_aux' = `r(mean)'
	
	sum transfer if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel O`counter_aux' = `r(mean)'
	
	sum transfer_pred_fgnv_lev_joint if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel P`counter_aux' = `r(mean)'
	
	sum tax_pred_fgnv_lev_joint if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel Q`counter_aux' = `r(mean)'
	
	sum snap_impute_val_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel R`counter_aux' = `r(mean)'
	
	sum housing_assist_impute_val_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel S`counter_aux' = `r(mean)'
	
	sum incwelfr_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel T`counter_aux' = `r(mean)'
	
	sum credits_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel U`counter_aux' = `r(mean)'
	    
	sum state_credits_hh if inc_tot_hh_groups == `i_group' [fw = asecwth]
	putexcel V`counter_aux' = `r(mean)'
	
	sum inc_lab_hh if inc_lab_hh_groups == `i_group' [fw = asecwth]
	putexcel W`counter_aux' = `r(mean)'
	
	sum tax_pred_fgnv_lev if inc_lab_hh_groups == `i_group' [fw = asecwth]
	putexcel X`counter_aux' = `r(mean)'
	
	sum tax if inc_lab_hh_groups == `i_group' [fw = asecwth]
	putexcel Y`counter_aux' = `r(mean)'

}

putexcel save

scalar drop _all